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Abstract 

The Boltzmann-Langevin dynamics of harmonic modes in nuclear matter is an- 
alyzed within linear-response theory, both with an elementary treatment and by 
utilizing the frequency-dependent response function. It is shown how the source 
terms agitating the modes can be obtained from the basic BL correlation kernel 
by a simple projection onto the associated dual basis states, which are propor- 
tional to the RPA amplitudes and can be expressed explicitly. The source terms 
for the correlated agitation of any two such modes can then be extracted directly, 
without consideration of the other modes. This facilitates the analysis of collec- 
tive modes in unstable matter and makes it possible to asses the accuracy of an 
approximate projection technique employed previously. 
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1 Introduction 



Microscopic transport models are necessary for the interpretation of heavy-ion 
colhsion experiments. Approaches based on the nuclear Boltzmann equation have 
been especially successful in accounting for a variety of observables, including 
the collective flow pattern and the production of mesons Models of this 

type seek to describe the one-body phase-space density f{r,p,t) of the nucleons 
in the system (and of any other hadron species present). The single-particle 
motion is governed by an effective Hamiltonian, h[f] = /2m + U{r), and the 
corresponding coUisionless evolution of / is governed by the Vlasov equation 0], 
the semiclassical analog of the Time-Dependent Hartree-Fock equation [Q. 

The advantage of the semiclassical description is that the constituent particles 
can be localized in phase-space and so it is possible to take account of their 
residual short-range interaction by means of a collision term, as was first done by 
Boltzmann for dilute classical gases. Due to the fermion nature of the nucleons, 
the two-body scattering processes are suppressed by the blocking factors / = 
1 — / expressing the availability of a state. This approach was first taken by 
Nordheim for the study of electrons in solids p, Q and was later adapted for 
nuclear collision scenarios |0, ||. This type of description is most often referred 
to as the Boltzmann- Uhling-Uhlenbeck {BUU) model 

In the BUU approach, only the average effect of the collisions are taken into 
account. While this suffices for many phenomena, it is inadequate for processes 
exhibiting bifurcations, such as nuclear fragmentation processes. The description 
has therefore been extended to include the fluctuating effect of the two-body 



collisions, leading to the Boltzmann-Langevin (BL) model [T^. The equation of 
motion for / is then of the form 

where I{r,p,t) represents the average effect of the collisions (the BUU term), 
while SI{r,p,t) denotes the fluctuating remainder (the Langevin term). 

The presence of the stochastic Langevin term 61 in the equation of motion 
for / causes continual branchings of the dynamical trajectories, thus enabling the 
system to explore the entire range of possible fates and, in particular, to exploit 
any instabilites encountered. The evolution of the corresponding distribution 
of histories, (f)[f]{t), can be described by a Fokker-Planck transport equation 



11, 12]. The associated transport coefficients, the drift coefficient V and the 



diffusion coefficient D, are given by 

V[f]{r,p,t) =^I{r,p,t) y , (2) 
2b[f] (r, p, p', t)5ir - r')5{t -t') = ^ 5/(r, p, t)5I{r\ p, t')* ^ , (3) 

where -< ■ >- denotes an average over an ensemble of systems that have all been 
prepared to have the same one-body density / at the time t. The transport coef- 
ficients are useful because they express how a given phase-space density f{r,p) 
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changes over a short period of time as a result of the collsions: the drift coefficient 
is the average rate of change, -< f >~, while the diffusion coefficient expresses the 
early growth rate of the correlation between the changes at the specified phase- 
space points, a{r,p;r',p'). Simple analytical approximations were recently de- 
rived for these key quantities 

The Langevin extension of the nuclear Boltzmann model is a significant ad- 
vance, since it makes it possible to address processes in which fiuctuations play 
a major role, such as nuclear multifragmentation caused by the irreversible de- 



velopment of unstable bulk modes [H. To gain insight into this key process, a 
recent study addressed the early evolution of nuclear matter in the spinodal zone 
of the phase diagram |T^. The system is then mechanically unstable and the 



density fluctuations generated by the Langevin term 61 may be amplified by the 
self-consistent effective field, leading towards catastrophic transformations of the 
system into an assembly of nuclear clusters. That work developed a convenient 
formal framework which provides instructive insight into the unstable dynamics 
and makes it possible to obtain quantitative results. 

It was shown that the collective modes are governed by simple feed-back 
equations of motion, in which the fiuctuations generated by the Langevin term 
are either magnified or suppressed by the adjusting effective field. The agitation 
rate is given in term of source terms V for which approximate expressions were 
derived, and they have been used to obtain quantitative results [1^. It is 



especially important to know accurately the agitation rates for the most unstable 
modes, since these will tend to become dominant, and in fact the final outcome 
may depend sensitively on their value. For example, an expanding system passing 
through the unstable spinodal zone may multifragment if sufficiently large density 
fiuctuations are produced during the finite time available, but if the fiuctuations 
remain too small the system may instead emerge relatively intact and then merely 
vaporize Moreover, if the fiuctuation amplitude is well understood, it may 
be possible to develop simplified stochastic mean-field approaches, where the 
complicated Langevin term is replaced by a simpler approximate term leading to 
the same agitation rates [jl8|, |19| . 



Therefore, in the present paper, we revisit this problem. We shall first address 
the situation by elementary means, thereby gaining instructive insight into the 
key problems, in particular the nature of the approximate treatment employed 
earlier [|l5l. Subsequently, we employ the frequency- dependent response function 
to obtain the same result in a more direct manner. As a main result, we find the 
explicit form of the dual basis states which makes it possible to obtain the collec- 
tive source terms directly. Furthermore, the accuracy of the inital approximate 
treatment given in ref. ITH] is illustrated. 



2 Elementary treatment 

We are considering the early evolution of nuclear matter that has been prepared 
with a uniform density distribution po and with a Fermi-Dirac momentum distri- 
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bution corresponding to a specified temperature T. It is therefore convenient to 
consider the Fourier transform of the phase-space density, 

/k(p,t)= J^e-^'^-'' Sfir,p,t) , (4) 



where 6f = /— /o is the deviation of f{t) from the initial density /o which depends 
only on the energy e. The system is confined within a volume Q, with periodic 
boundary conditions imposed, so that the wave number k takes on discrete values. 
As long as we remain within the regime of linear response, the dynamical problem 
decouples with respect to wave number, and so we may limit our considerations 
to a single value of k. 

The equation of motion for these Fourier coefficients follows readily from (|lD, 

— /k + ik ■ v/k - ■ v^p\, = dik , (5) 

Here dUk/dp represents the appropriate Fourier component of the derivative 
of the effective field with respect to the density and 5/k(p, t) is the Fourier 
component of the Langevin term in (|l|), defined in analogy with eq. (^). Fur- 
thermore, pk(i) = J2p fiiiP^t) is the Fourier transform of the matter density 
p(^) = Hpf{fiP)i where X]p denotes the summation over all the momentum 
states and is equivalent to the integral J dp in the limit of large systems 
having a continuum of momentum states.^ Finally, we have ignored the average 
collision term J, since its effect is relatively small, consisting primarily of a small 
reduction in the growth times for the unstable modes |^D|, pT| . 

The Vlasov equation (the left-hand part of (^) has a complete set of eigen- 
functions, 

/k(p,t) = /k(p)e-^'^''*, (6) 

on which any time-dependent solution to may be expanded. These eigenfunc- 
tions satisfy the equation 

/I \pu/ \ 9Uk dfo 

{k-v-uj^)f^{p) = -^k-v— p^ , (7) 

and so they have the form 

dUk kv dfo 



KiP) 



dp k - V — tUu de 



where a specific normalization has been chosen, p'^ = J2p fkip) = 1- 
The eigenfrequencies Ui, are determined by the dispersion relation, 

dUk^ kv dfo f ^ /n^ 

dp p k - V — uOy de 



^We prefer to employ the summation sign, since the spectral structure is easier to visualize 
when the problem is discretized by means of a lattice in momentum space. 
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In the last expression, we have employed the usual Landau parameter Fq and the 
free response function 

which is the average of k ■ v/(k ■ v — cu) over momentum space, calculated with 
the weight function dfo/de = —fofo/T emphasizing the Fermi surface. It is easy 
to see that the solutions come in pairs of opposite signs (which we may assign 
opposite values of the index u: uj^u = —uju). 

In the absence of the damping term /, the dispersion relation (P) has ex- 
clusively real roots as long as the system is prepared outside the spinodal zone, 
i.e. when Fq > —1. However, inside the spinodal zone one pair of eigenvalues is 
purely imaginary, ujk = ^i/tk, and is associated with a pair of collective modes, 
one exponentially growing and the other one exponentially decaying, with the 
same characteristic time constant t^. The amplified collective mode will quickly 
become the dominant one and consequently this is the mode of primary concern. 
However, our present developments are quite general and apply to all the modes, 
and they are not limited to the spinodal zone. 

It is important to recognize that although the eigenfunctions {f^} form a 
complete set, they are not orthogonal. We therefore introduce the matrix or as 



the inverse of the overlap matrix [15 



iok'r' =< E/k(p)Vk^'(p) • (11) 

p 

Becuase of the rotational invariance, the overlap matrix elements do not depend 
on the direction of the wave vector k but only on its magnitude k. Any solution to 
the equation of motion (||), /k(p, t), has a unique expansion on the eigenfunctions, 

/k(p,t) = E^kW/k(p) ' (12) 

V 

where the expansion coefficients are given by 

K = Y.o7' < K\h> = < €\h> ■ (13) 

We have here introduced the functions ^^(p) = J2u' o^"' fk (p) "which form the 
dual basis, relative to the eigenfunctions {f^{p)}, since 

< 9kl/k' > = , (14) 



as is readily shown by using the definition ([11]) of the overlap matrix. 

Inserting the expansion (^) into the equation of motion (^ and projecting 
onto the dual basis we obtain the following equation for the expansion amplitudes, 

j^A'^ = -zu;,A'^{t)+ < gl^|<5/k(t) > . (15) 
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This form is easy to understand: the noise term SI acts as a source term that 
continually produces stochastic changes of the amplitude A"^, while the Vlasov 
equation propagates the amplitude in the associated effective field. The above 
equation of motion can readily be solved formally, 

AUt) = fdt' e^-^(*'-*) < ql\5Ut') > + AUO) . (16) 
Jo 

If we start from uniform matter, the initial amplitudes vanish, ^^(0) = 0, and 
since the noise term vanishes on the average, -< 61 >~= 0, the ensemble averages 
of the amplitudes remain zero, -< A'^{t) >~= 0. However, their correlations have 
a non-trivial evolution, 

-< A'C(t)AC(ty y = 5kk' 2Vf t dt' e'(— :')(*'-*) . (17) 

Jo 

We have here used the fact that the Langevin term is local in space and time, so 
the problem decouples with respect to the wave number k, 

< 6h{p,t)6hip',ty >- = 2D{p,p') 6{t - t') 4k' . (18) 

Furthermore, we have introduced the source term 

Vf = < q\^\D\q'C > = E <lM* D{p,p') q^p') . (19) 

pp' 



We note that the above result (|T^ is equivalent to the following equation of 
motion for the correlation coefficients a^^ (t) =^ A'^{t)A'(^ (t)* y, 

j^afit) = 2Vf - - a.:,)ar'(t) , (20) 



which is the result derived in ref. [|15|. In order to better make contact with that 
work, we note that the expression (|T9|) for the source term can be expressed in 
terms of the basis formed by the eigenf unctions, 

^^r' = < €\D\qi > = E«r < f^m( > of ■ m) 

In ref. the focus was on the evolution of the two collective modes in the 
spinodal zone and only those modes were retained. So the full matrix o^'^ was 
effectively replaced by the 2x2 matrix 5^'^ involving only the two collective 
modes. This amounts to approximating the exact dual state ^{^(p) by gk(p) = 
/k (p)) where the sum extends only over the two collective modes. The 
corresponding approximate expression for the collective source term is then 

w' = < iii\D\€ > = E^r < fmf^' > of , (22) 



as given in ref. []T5|. This analysis exhibits the relationship between the two 
projection methods, the initial one (^) using the approximate dual basis func- 
tions and the improved one ( PT| ) employing the exact dual basis. In section ^ we 
illustrate this further by numerical comparisons. 
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3 Response function approach 



It is possible to treat the problem in an alternative manner, by direct application 
of the response-function technique as is briefly summarized below. 

Considering a general solution to (H), fk,{p,t), we first perform a one-sided 
Fourier transform with respect to time, 

Mp,uj)= / dt e'-' f^{p,t) , (23) 
Jo 

leading to the equation 

i(k • V - t^)/k(p, uj) - ^-^k ■ v-^p]^{uj) = 51]^{p, u) , (24) 

where SIy_(p,uj) is the corresponding transform of 6I\f_{p,t). 
It is common to introduce the susceptibility 

e{uj) = 1 + FoxH = 1-^2^ ] TT- , 25 

op p K - V — CO oe 

which is simply related to the response function x{^) introduced in ([To|) . The 
solution of the above equation (|2^) can then be written as a simple projection of 
the noise, 

p e{uj) p K - V — 00 

The Fourier component of the density can subsequently be obtained by applying 
the inverse transform, 

f u,i / N r , , f dio e^'^(*'-*) 5Jk(p,t') 

P^(t) = L.^ P^(-^) = L dt' I — -—— E r^T^ • (27) 



c 271 Jo Jc 2ni e^oo) p k ■ v — oo 

Here the integral is to be carried out along a path C in the complex plane that 
passes above all the poles of the integrand, from the far left to the far right. When 
t' > t the Lo integral can be carried out by completing the contour around the 
upper halfplane where the integrand vanishes rapidly. Since the completed con- 
tour encloses no poles, the uo integral vanishes and so the t' integration effectively 
extends only up to t' = t. To evaluate the u integral for t' < t, we complete the 
contour around the lower halfplane, which adds no contribution, and then apply 
the Residue Theorem. The poles are determined by the condition e{u^) = 0, 
which is recognized as the dispersion relation (P) and so they occur at the eigen- 
frequencies, as expected. Furthermore, the residue of the susceptibiltiy is 

\ / LiJ=UJv P ^ 
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where we have introduced the auxihary function Q'^{p) = l/(k • v — cj*), which 
is recognized as the usual RPA amphtude. We may then write the result (|27|) on 
a compact form, 

P.{t) = f dt' E e-(*'-') = E Kit) . (29) 

The last relation arises by recalling the expansion (|1^) of the solution on the 
eigenfunctions, and the fact that we use a normalization such that p'^ is unity. 

Because of the stochastic character of the evolution, we need to consider the 
entire ensemble of possible dynamical histories. The ensemble average of the 
Langevin noise term vanishes, -< 5/k(p, t) >~= 0, and so -< A'^(t) >~= 0, as in sect. 
1^. The corresponding correlation function is given by 

The results (|^) and ( pOf ) are identical to what was derived in sect. H, since the 
auxiliary function Q\.{p) is proportional to the normalized dual basis function 

This key relation can be verified by using the dispersion relation (j^) and the fact 
that the eigenvalues come in opposite pairs, as is demonstrated in the Appendix. 

The above derivation shows that the Boltzmann-Langevin dynamics in nuclear 
matter can be treated by suitable adaptation of the standard response-function 
approach pO| , p^ . The availability of an explicit form of the dual basis enables 
us to determine the source term directly for each pair of modes v and z/', 
without the need for invoking the entire overlap matrix o'^^y . This is particularly 
useful inside the spinodal zone, since the dynamics quickly becomes dominated 
by the amplified collective mode, so the there is little need for considering the 
others. Accordingly, the dual projection method was employed in a recent study 
of the effect of memory time on the agitation of unstable modes in nuclear mattter 

m. 



4 Discussion 

We have shown how the Boltzmann-Langevin dynamics of nuclear matter can 
be treated within linear response theory. This treatment is valid as long as the 
induced deviations remain small and so, for matter in the unstable phase region, 
it can be used to understand the onset of the spinodal decomposition. 

The time evolution of the density undulations can be separated into two qual- 
itatively different stages. The characteristic time separating the two dynamical 
regimes is given by the amplification time tk. Whereas all the modes are initially 
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agitated in a rather democratic manner, in accordance with the miscoscopic diffu- 
sion coefficient D(j)^,P2)^ only the amphfied collective mode will ultimately dom- 
inate (since the suppressed collective mode will saturate and the non-collective 
ones have no exponential development). Therefore, if one is interested in the 
long-term behavior, i.e. the appearance of the system after a time longer than 
the characteristic time, then it suffices to retain only the source term for the 
amplified collective mode. That is why we have focussed on the dynamics of the 
collective modes. The agitation rate for the amplified collective mode is most eas- 
ily obtained by projecting the diffusion coefficient D{pi,P2) onto the dual basis 
state for that mode, q^{p), as indicated in eq. (pOl). 

However, it should be recognized that for times short in comparison with the 
amplification time, t < tk, the relative magnitude of the contributions from other 
modes is not yet negligible. This is most clearly brought out by the fact that the 
non-collective contributions cancel exactly the collective ones to leading order in 
time, so that ak ~ t"^ at first, whereas the expression (^) yields a linear initial 
growth. This general property follows formally from particle-number conserva- 
tion which dictates that the integral of the microscopic diffusion coefficient with 
respect to momentum vanish. The feature can easily be understood from the 
fact that the collision term is local in space and time. It then produces merely 
a rearrangement of the local occupation in momentum space, without affecting 
directly the matter density, and fluctuations in the matter density only appear 
subsequently as the generated two-particle-two-hole excitons are propagated in 
the effective field. Consequently the effect of the Langevin term on the matter 
density is only second order in time. Contrary to this, the density fluctuation 
arising from the agitation of the amplified mode alone is initially linear in time. 

The key quantities, such as the collective source terms, can be evaluated 



either by direct numerical calculation [15, 16] or by employing the analytical 



approximation derived for the diffusion coeffision in ref. |jT3[ and making use of 
the dispersion relation. For our present illustrative purposes, we follow the latter 
approach, which is discussed in detail in ref. and employs a realistic two-body 
interaction for the calculation of the Landau parameter Fq. 

Let us first recall that the spinodal zone is situated within a parabola-like 
boundary that extends from zero density to about two thirds of the saturation 
density for vanishing temperature and reaches upwards to the critical temperature 
of about 16 MeV for symmetric and uncharged nuclear matter. This standard 
spinodal boundary pertains to undulations with a vanishing wave number, k = 0, 
and it shrinks as the wave number is increased, and finally disappear altogether 
at a point located on the T = axis near one third of the saturation density. The 
amplification time tk{p, T) is infinitely large at the boundary and becomes shorter 
as the phase point considered is moved closer towards that limiting point. The 
fastest amplification time is about 0.7 • 10~^^ s but typical values are somewhat 
larger, due to the finite temperature. 

Figure |l] shows the collective source term T)^^ obtained with the model pre- 
sented in ref. PT|. We note that this quantity diverges at the spinodal boundary 
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where Fq — > —1, which may serve as a reminder of the hmitations mentioned 
above. Near the boundary the characteristic time tends to infinity, ^ oo, and 
so we remain in the early dynamical regime where the other modes are signifi- 
cant; indeed, they will conspire to cancel exactly the divergent behavior of the 
amplified mode. 

The original orthogonal projection introduced in ref. effectively employs 
an approximate dual basis function obtained by decoupling the 2x2 collective 
part of the overlap matrix from the rest. The result is shown in fig. |l| by the 
dashed curves and it is seen to behave in a manner very similar to that obtained 
by projecting onto the exact dual basis function (solid curves). Indeed, it exhibits 
a similar divergent behavior near the spinodal boundary. However, it is typically 
smaller by about 30% for the most rapidly amplified modes, and so the use of 
the exact dual basis represents a significant quantitative improvement. 

In order to compare the two methods with regard to their predictions for the 
density fiuctuations, we consider the Fourier component of the density variance. 



where the last expression has been obtained by retaining only the two collective 
modes, and it has been used that V'^^ = Vj^' and = T>j^^. The resulting 
variance is displayed in fig. |^ for a sequence of times t. The approximate projec- 
tion method involves the inversion of the 2x2 collective overlap matrix. As the 
spinodal boundary is approached, the collective frequency tends to zero and the 
two collective modes become identical. The mixed source term then also diverges 
and, as it happens, V^^ + V^^ ~ t^-^ as tjt oo. So the resulting density 
variance ak tends to zero. Such a regular behavior does not arise when the exact 
projection is employed, because the mixed source term remains regular and so 
cannot cancel the divergence of the diagonal term. The resulting density variance 
then diverges at the spinodal boundary, emphasizing the fact that the method is 
limited to scenarios well within the spinodal zone where the characteristic time 
is reasonably short. 

In fig. 1^ we show the density variance ak as a function of the wave number 
k = 27r/A and at the same times, for the typical temperature T = 4 MeV and 
for a density near which the fastest amplification occurs, p = 0.3po- It is seen 
how those modes that have the shortest amplification time grow progressively 
dominant, so that the density fiuctuations emerge with a rather narrow Fourier 
spectrum centered around the fastest-growing mode associated with the specified 
density and temperature. It should be noted that the exponential growth limits 
the length of time over which the linear-response treatment remains valid. From 
the results shown in fig. ^ it can readily be estimated that the predicted average 
magnitude of the density fiuctuation equals the average density for t ~ 3.5-10"^^ s, 
so the results are not expected to be meaningful beyond t 3 ■ 10~^^ s. A 
comparison between the results of the two different projection methods leads to 




-< pIp^ y ^ 2V++tu sinh(-) + AV^-t , 



(32) 
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the same conclusion as reached above: use of the approximate dual basis states 
reduces the density variances by about 30%, except near the spinodal boundary 
where the characteristic time diverges and it becomes insufficient to consider only 
the collective modes. 

While we expect that the projection method employing the exact dual basis 
states is quite reliable for times exceeding the respective amplification time, t > 
tk, it still needs to be determined how accurate it is at earlier times. In particular, 
it would be of interest to improve the projection method so that it is applicable 
also near the spinodal boundary. This task is of practical importance, since the 
systems prepared in heavy-ion collisions are initially situated outside the spinodal 
region and so must cross the boundary to become unstable. 
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A The explicit form of the dual basis 

In this Appendix we prove that the auxiliary functions Qk(p) ~ 1/ (k-v — cj*) are 
indeed proportional to the dual basis q^{p) characterized by < g^l/k >= ^fi^'- 
We ffist recall that the eigenvalues come in pairs of opposite sign and that 
we label them accordingly, u^i, = —u^. Let us start with the most frequent case 
when the two eigenvalues differ, cUi, 7^ u^/, 

<wr>-E ,. ' .... t (33) 



p 



]i ■ V — uji,\i ■ V — ujyi de 
k v dfo V- k ■ V dfo 



uji, — ujyi I p k - V — ujy at pk-v — ujyi ot I 

This result is obtained after a few of elementary manipulations, followed by appli- 
cation of the dispersion relation (^ to each of the two terms, which then cancel. 
This leaves only the possiblity the the two eigenvalues are equal, z/ = z/', in which 
case we have 

< «|/- >= Y. , ' ^ = E , ^ . (35) 

^kuk ^k-v-uj^k-v-uj^ de ^ {{k ■ vf - ujlf de ^ ^ ' 

Therefore, after a suitable renormalization, = A/kQ^^ ^^e required orthonor- 
mahty relation follows, < gkl/k >= ^w'-. ^-^^ {O'k} indeed forms the dual basis. 
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Figure 1: Comparison of source terms. 
The source term T>^^ for the amphfied collective mode having a wave length A = 
27r/fc, in nuclear matter prepared with a uniform density p and with a specified 
temperature T. The solid curves are obtained by projecting onto the exact dual 
basis state, eq. (0), while the dashed curves show the corresponding results 
obtained with the approximate dual basis state, eq. The calculations have 



been done with the approximate formulas developed in ref. [21 . The upper panel 



considers the typical temperature T = 4 MeV and illustrates the dependence 
on the wave length A, while the lower panel keeps the wave length fixed at A = 
27r//c = 8 fm, near which the most rapid amplification occurs, and illustrates the 
temperature dependence. 



Figure 2: Density fluctuations. 
The variance of the density fluctuations associated with the collective modes, ak, 
after a given time t = 1, 2, 3, 4 ■ 10"^^ s has elapsed, as a function of the average 
density and for the temperature T = 4 MeV and a wave length of A = 8 fm. 
Notation and model are as in fig. |l]. 



Figure 3: Dependence on wave number. 
The variance of the density fluctuations associated with the collective modes, ak, 
after a given time t = 1,2,3,4 ■ 10~^^ s has elapsed, as a function of the wave 
number k = 2n/X, for the typical temperature T = 4 MeV and for the density 
p = 0.3po- Otherwise similar to fig. ^ 
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Figure 1b 
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